
#include "main.h"
#include "../mds/engine.h"
#include "../mds/constants.h"
#include "../mds/hasten/cells-dyn.h"
#include <math.h>
#include "../mds/io/xyz.h"
#include "../mds/io/mds.h"
#include "../mds/sys/string.h"
#include "../mds/density.h"
#include "../mds/sys/pid.h"

/*******************************************************************************
 *
 * 
 *
 ******************************************************************************/
#define MD_ERROR exit_value = -1; goto MD_END;

#define on_option_initial     on_option_id
#define on_option_restart     on_option_id
#define on_option_output      on_option_id
#define on_option_box         on_option_id
#define on_option_mass        on_option_id
#define on_option_T           on_option_id
#define on_option_benchmark   on_option_id
#define on_option_np          on_option_id 

/*******************************************************************************
 *
 * saving
 *
 ******************************************************************************/
static
int save_energy( const Real t, const mds_frame *frame, int append)
{
	FILE *fp = fopen( "energy.txt", append ? "ab" : "wb" );
	Real T=0, PV=0;
	if( !fp )
		return -1;
	mds_frame_NVE_T_and_PV(frame, &T, &PV);
	fprintf( fp, "%.6e %.6e %.6e %.6e %.6e %.6e %.6e\n", t, frame->V, frame->K, frame->V+frame->K, frame->W, T, PV );
	fclose(fp);
	return 0;
}

/*******************************************************************************
 *
 * main
 *
 ******************************************************************************/

int main( int argc, char *argv[] )
{
	int            exit_value = 0;
	mds_step_type  step_type  = MDS_STEP_VERLET;
	mds_engine     engine     = MDS_ENGINE_INIT;
	mds_options    options    = MDS_OPTIONS_INIT;
	Real           t_ini      = 0;
	Real           t_run      = 0;
	size_t         i_run      = 0;
	Real           t_max      = 0;
	size_t         every      = 0;
	const char    *initial    = NULL;  /*!< load initial frame: will save first frame    */
	const char    *restart    = NULL;  /*!< load restart frame: won't save first frame   */
	const char    *output     = NULL;  /*!< output file name: can be restart             */
	int            genconf    = 0;     /*!< if no initial or restart                     */
	const char    *np_str     = NULL;  /*!< #particles if no initial|restart             */
	size_t         np         = 0;     /*!<  #particles if no initial|restart            */
	const char    *box_str    = NULL;  /*!< not null => override configuration           */
	const char    *mass_str   = NULL;  /*!< not null => 1..mass if no initial|restart    */
	const char    *T_str      = NULL;  /*!< not null => override T if no initial|restart */
	const char    *benchmark  = NULL;  /*!< not null => benchmark */
	const char    *benchfile  = "benchmark.txt";
	
	
	mds_verlet_data verlet_data;
	mds_beeman_data beeman_data;
	
	mds_ljc_data ljc_data;
	Real         epsilon = 1;
	Real         sigma   = 1;
	Real         T       = 0;
	
	mds_vector *position = NULL;
	mds_vector *velocity = NULL;
	
	mds_dynamics     whole_dynamics;
	mds_dynamics     crude_dynamics;
	mds_dyn_type     dyn_type = MDS_DYN_CELL;
	mds_run_type     run_type = MDS_RUN_NVE;
	int              rescaleV = 0;
	Real             K0;
	mds_wtime        chrono;
	mds_wtime        watch;
	double           crunch = 0;
	register size_t i;
	
	/*==========================================================================
	 *
	 *  Initializing
	 *
	 *========================================================================*/
	mds_run_initialize();
	MDS_R = 1;
	engine.frame.pbc_for_f = mds_pbc_xyz;
	engine.frame.pbc_for_r = mds_pbc_xyz;
	mds_wtime_init( &chrono );
	mds_wtime_init( &watch  );
	
	/*==========================================================================
	 *
	 *  Options Loading
	 *
	 *========================================================================*/
	fprintf( stderr, "-----------------------------------------------------\n");
	fprintf( stderr, "-- Loading Options \n");
	fprintf( stderr, "-----------------------------------------------------\n");
	MDS_OPTIONS_DECL()
	{
		MDS_OPTION_NORMAL(f,         options,          "paramaters file"),
		MDS_OPTION_FORCED(dt_sim,    engine.dt_sim,    "simulation time step"),
		MDS_OPTION_FORCED(dt_out,    engine.dt_out,    "ouput time step"     ),
		MDS_OPTION_FORCED(t_max,     t_max,            "simulation time max" ),
		MDS_OPTION_FORCED(step,      step_type,        "simulation step : verlet|beeman "),
		MDS_OPTION_FORCED(cutoff,    engine.cutoff,    "potential cut off length"),
		MDS_OPTION_NORMAL(dyn,       dyn_type,         "dynamic type: full|cell (default=cell)"),
		MDS_OPTION_FORCED(cutout,    engine.cutout,    "cells  cut out length" ),
		MDS_OPTION_NORMAL(initial,   initial,          "load initial frame from file (start new simulation)"),
		MDS_OPTION_NORMAL(restart,   restart,          "load restart frame from file (continue simulation)"),
		MDS_OPTION_NORMAL(output,    output,           "output file, can be restart" ),
		MDS_OPTION_NORMAL(np,        np_str,           "#particles, if no initial|restart file"),
		MDS_OPTION_NORMAL(box,       box_str,          "override if initial|restart, mandatory if no initial|restart"),
		MDS_OPTION_NORMAL(mass,      mass_str,         "random masses in [1;mass] if no initial|restart (default=1)"),
		MDS_OPTION_NORMAL(T,         T_str,            "override f initial|restart, mandatory if no initial|restart"),
		MDS_OPTION_FORCED(run,       run_type,         "run NVE|NVT"),
		MDS_OPTION_NORMAL(benchmark, benchmark,        "enter benchmark mode / #bench times")
	}
	MDS_OPTIONS_DONE();
	
	if( argc <= 1 )
	{
		MDS_OPTIONS_HELP(stderr);
		exit_value = 1;
		goto MD_END;
	}
	
	if( mds_options_parse( &options, argc, argv) < 0 )
	{
		mds_perror( stderr, "mds_options_parse()");
		MD_ERROR;
	}
	
	MDS_OPTIONS_DISPATCH(options,MD_ERROR);
	
	/*==========================================================================
	 *
	 *  Options Checking
	 *
	 *========================================================================*/
	fprintf( stderr, "-----------------------------------------------------\n");
	fprintf( stderr, "-- Checking Options \n");
	fprintf( stderr, "-----------------------------------------------------\n");
	
	/*--------------------------------------------------------------------------
	 *-- time steps
	 *------------------------------------------------------------------------*/
	if( engine.dt_out <= engine.dt_sim ) engine.dt_out = engine.dt_sim;
	
	fprintf( stderr, "mdc> dt_sim=%g | dt_out=%g\n", engine.dt_sim, engine.dt_out );
	every = (size_t)(floor(engine.dt_out/engine.dt_sim));
	if( every <= 0 ) every = 1;
	engine.dt_out = every * engine.dt_sim;
	fprintf( stderr, "mdc> saving every %u steps, dt_out=%g\n", (unsigned)every, engine.dt_out );
	
	/*--------------------------------------------------------------------------
	 *-- I/O
	 *------------------------------------------------------------------------*/
	if(output) 
		fprintf( stderr, "mdc> output to '%s'\n", output );
	else 
		fprintf( stderr, "mdc> NO OUTPUT !!!\n");
	
	if( restart && initial )
	{
		exit_value = -1;
		fprintf( stderr, "** can't both restart and start from initial!\n");
		goto MD_END;
	}
	
	if( initial )
		fprintf( stderr, "mds> loading initial frame from '%s'\n", initial );
	
	if( NULL == initial && NULL == restart )
	{
		genconf = 1;
		fprintf( stderr, "mds> create a new configuration\n");
	}
	
	/*--------------------------------------------------------------------------
	 *-- Choose Step
	 *------------------------------------------------------------------------*/
	switch( step_type )
	{
		case MDS_STEP_VERLET:
			fprintf( stderr, "mdc> VERLET step\n" );
			mds_verlet_step_load( & engine.step, &verlet_data );
			break;
			
		case MDS_STEP_BEEMAN:
			fprintf( stderr, "mdc> BEEMAN step\n" );
			mds_beeman_step_load( & engine.step, &beeman_data );
			break;
	}
	
	assert( engine.step.num_sets > 0 );
	
	/*--------------------------------------------------------------------------
	 *-- cutoff, epsilon, sigma...
	 *------------------------------------------------------------------------*/
	fprintf( stderr, "mdc> initializing potential\n" );
	if( engine.cutout <= engine.cutoff )
	{
		fprintf( stderr, "** cutout=%g <= cutoff=%g\n", engine.cutout, engine.cutoff);
		MD_ERROR;
	}
	
	mds_ljc_init( &ljc_data, epsilon, sigma, engine.cutoff);
#if 0
	{
		FILE *pot = fopen("vf.txt", "wt");
		if( pot )
		{
			Real xx;
			for( xx = 0.95*sigma; xx <= 2 * engine.cutoff; xx += 0.001 )
			{
				mds_vector R    = { xx, 0, 0 };
				mds_vector F    = { 0, 0,  0 };
				const Real V    = mds_ljc( &R, &F, &ljc_data );
				const Real V_LJ = mds_lj( &R, &F, &ljc_data );
				fprintf( pot, "%g %e %e %e\n", xx, V, F.x, V_LJ );
			}
			fclose( pot );
		}
	}
#endif
	
	
	/*==========================================================================
	 *
	 *  Generate Configuration Or Load it
	 *
	 *========================================================================*/
	if( genconf )
	{
		/*----------------------------------------------------------------------
		 *
		 *  generating configuration
		 *
		 *--------------------------------------------------------------------*/
		mds_uvector cr;
		Real        mass_max = 1;
		
		fprintf( stderr, "-------------------------------------------------\n");
		fprintf( stderr, "-- Generating Configuration \n");
		fprintf( stderr, "-------------------------------------------------\n");
		/*----------------------------------------------------------------------
		 *  we need particles
		 *--------------------------------------------------------------------*/
		if( !np_str || (np=strtol(np_str,NULL,10)) <= 0 )
		{
			fprintf( stderr, "** no particle to generate configuration!\n");
			MD_ERROR;
		}
		fprintf( stderr, "mdc> #particles=%u\n", (unsigned)np );
		
		/*----------------------------------------------------------------------
		 *  we need a box
		 *--------------------------------------------------------------------*/
		if( ! box_str )
		{
			fprintf( stderr, "** no box to generate configuration !\n");
			MD_ERROR;
		}
		
		if( mds_box_parse( & engine.frame.box, box_str ) < 0 )
		{
			fprintf( stderr, "** invalid box format '%s'!\n", box_str);
			MD_ERROR;
		}
		
		fprintf( stderr, "mds> generating box = [ %g %g %g ]\n", engine.frame.box.x, engine.frame.box.y, engine.frame.box.z );
		
		/*----------------------------------------------------------------------
		 *  Check for mass
		 *--------------------------------------------------------------------*/
		if( mass_str )
		{
			mass_max = strtod( mass_str, NULL);
			if( mass_max < 1 ) 
			{
				fprintf( stderr, "** invalid mass = %g !\n", mass_max);
				MD_ERROR;
			}
		}
		
		/*----------------------------------------------------------------------
		 *  Check for temperature
		 *--------------------------------------------------------------------*/
		if( !T_str )
		{
			fprintf( stderr, "** no Temperature to generate configuration!\n");
			MD_ERROR;
		}
		
		T = strtod( T_str, NULL);
		if( T < 0 )
		{
			fprintf( stderr, "** invalid Temperature!\n");
			MD_ERROR;
		}
		fprintf( stderr, "mdc> initial Temperature= %g\n", T );
		
		/*----------------------------------------------------------------------
		 *  now we count them on a crystal
		 *--------------------------------------------------------------------*/
		np = mds_crystal_close_packed_count(&cr, np, & engine.frame.box );
		fprintf( stderr, "mdc> #particles=%u on this fcc crystal\n", (unsigned)np );
		
		/*----------------------------------------------------------------------
		 *  Ready for the engine
		 *--------------------------------------------------------------------*/
		if( mds_engine_load( &engine, np) < 0 )
		{
			fprintf( stderr, "** can't load engine\n");
			MD_ERROR;
		}
		
		
		/*----------------------------------------------------------------------
		 *  Prepare mass/inv_mass and velocities
		 *--------------------------------------------------------------------*/
		position = engine.frame.sets[ MDS_POSITION ];
		velocity = engine.frame.sets[ MDS_VELOCITY ];
		for( i=0; i < np; ++i )
		{
			mds_vector *r    = &position[i];
			mds_vector *v    = &velocity[i];
			const Real  mass = 1 + (mass_max-1) * mds_alea();
			
			v->info.mass     = mass;
			r->info.inv_mass = ((Real)1) / mass;
		}
		
		/*----------------------------------------------------------------------
		 * put particles at their place
		 *--------------------------------------------------------------------*/
		mds_crystal_build_fcc(position, np, &cr, & engine.frame.box, 0);
		
		
		/*----------------------------------------------------------------------
		 *  prepare initial velocities
		 *--------------------------------------------------------------------*/
		mds_frame_initialize_velocities( &engine.frame, T );
		
		
		/*----------------------------------------------------------------------
		 *
		 *  Done generating configuration 
		 *
		 *--------------------------------------------------------------------*/		
	}
	else
	{
		/*----------------------------------------------------------------------
		 *
		 * ENTER: Loading configuration
		 *
		 *--------------------------------------------------------------------*/
		const char *filename = restart ? restart : initial;
		mds_frame   target   = MDS_FRAME_INIT;
		target.pbc_for_r     = mds_pbc_xyz;
		fprintf( stderr, "-------------------------------------------------\n");
		fprintf( stderr, "-- Loading Configuration \n");
		fprintf( stderr, "-------------------------------------------------\n");
		
		
		/*******************************/
		/** create a temporary frame **/
		/*****************************/
		assert( filename != NULL );
		fprintf( stderr, "mdc> loading from...'%s'\n", filename );
		if( mds_file_load( filename, &target,& t_ini ) < 0 )
		{
			mds_perror( stderr, "mds_file_load(%s)", filename);
			MD_ERROR;
		}
		np  = target.num_particles;
		fprintf( stderr, "mdc> loaded %u particles\n", (unsigned) np );
		
		fprintf( stderr, "-------------------------------------------------\n");
		fprintf( stderr, "-- Adjusting configuration \n");
		fprintf( stderr, "-------------------------------------------------\n");
		
		/*****************************/
		/** can't change particles **/
		/***************************/
		if( np_str )
		{
			fprintf( stderr, "** can't set #particles while reloading file!\n");
			MD_ERROR;
		}
		
		/********************************/
		/** can change the target box **/
		/******************************/
		if( box_str )
		{
			mds_vector usr_box = MDS_VECTOR_INIT;
			if( mds_box_parse( &usr_box, box_str) < 0 )
			{
				mds_frame_quit( &target );
				fprintf( stderr, "** invalid box '%s'\n", box_str );
				MD_ERROR;
			}
			fprintf( stderr, 
					"mdc> rescaling from [%g %g %g] to [%g %g %g]\n", 
					target.box.x, target.box.y, target.box.z,
					usr_box.x, usr_box.y, usr_box.z);
			mds_frame_rescale( &target, & usr_box );
			fprintf( stderr, "mdc> target box is now [ %g %g %g]\n", target.box.x, target.box.y, target.box.z );
		}
		
		/***********************************************/
		/** make the engine from the temporary frame **/
		/*********************************************/
		mds_vector_cpy( &engine.frame.box, & target.box );
		fprintf(stderr, "mdc> creating engine\n");
		if( mds_engine_load( &engine, target.num_particles ) < 0 )
		{
			mds_frame_quit( &target );
			MD_ERROR;
		}
		position =  engine.frame.sets[MDS_POSITION];
		velocity =  engine.frame.sets[MDS_VELOCITY];
		memcpy( position, target.sets[MDS_POSITION], np * sizeof(mds_vector) );
		memcpy( velocity, target.sets[MDS_VELOCITY], np * sizeof(mds_vector) );
		engine.frame.K = target.K;
		engine.frame.V = target.V;
		engine.frame.W = target.W;
		
		/*********************************/
		/** delete the temporary frame **/
		/*******************************/
		mds_frame_quit( &target );
		
		
		/*********************************/
		/** check frame                **/
		/*******************************/
		fprintf( stderr, "mdc> %u particles in [ %g %g %g ]\n", (unsigned)engine.frame.num_particles, 
				engine.frame.box.x, engine.frame.box.y, engine.frame.box.z);
		
		T = mds_frame_T( &engine.frame );
		fprintf( stderr, "mds> << T= %g >>\n", T);
		
		if( T_str )
		{
			T = strtod( T_str, NULL );
			if( T < 0 )
			{
				fprintf( stderr, "** invalid temperature !!\n");
				MD_ERROR;
			}
			K0 = ( MDS_F(1.5) * MDS_R * T) * np;
			fprintf( stderr, "mdc> << Scaling to T= %g, K0=%g >>\n", T, K0 );
			engine.frame.K = mds_velocity_rescale( velocity, np, (Real)sqrt(K0/engine.frame.K) );
			fprintf( stderr, "mdc> starting K=%g\n", mds_velocity_kinetic_energy(velocity,np) );
		}
		
		/*************/
		/** timing **/
		/***********/
		if( initial ) t_ini = 0;
		
		
		/*----------------------------------------------------------------------
		 *
		 * LEAVE: Loading configuration
		 *
		 *--------------------------------------------------------------------*/
	}
	
	assert( np       == engine.frame.num_particles );
	assert( position == engine.frame.sets[MDS_POSITION] );
	assert( velocity == engine.frame.sets[MDS_VELOCITY] );
	
	/*==========================================================================
	 *
	 *  Prepare Potential and Dynamics
	 *
	 *========================================================================*/
	fprintf( stderr, "-----------------------------------------------------\n");
	fprintf( stderr, "-- Adjusting potential and dynamics\n");
	fprintf( stderr, "-----------------------------------------------------\n");
	
	switch( dyn_type )
	{
		case MDS_DYN_FULL:
			whole_dynamics.gather = mds_dynamics_full;
			whole_dynamics.params = NULL;
			
			crude_dynamics.gather = mds_dynamics_full;
			crude_dynamics.params = NULL;
			
			benchfile = "benchmark-full.txt";
			fprintf( stderr, "mdc> using FULL dynamics\n");
			break;
			
		case MDS_DYN_CELL:
			whole_dynamics.gather = mds_cells_dynamics_whole;
			whole_dynamics.params =  & engine.cells;
			
			crude_dynamics.gather = mds_cells_dynamics_whole;
			crude_dynamics.params = & engine.cells;
			
			fprintf( stderr, "mdc> using CELL dynamics\n");
			benchfile = "benchmark-cell.txt";
			break;
	}
	
	K0 = engine.frame.K;
	switch( run_type )
	{
		case MDS_RUN_NVE:
			fprintf( stderr, "mdc> running in NVE\n");
			break;
			
		case MDS_RUN_NVT:
			fprintf( stderr, "mds> running in NVT\n");
			rescaleV = 1;
			break;
	}
	
	
	
	/*==========================================================================
	 *
	 *  preparing to run
	 *
	 *========================================================================*/
	fprintf( stderr, "-----------------------------------------------------\n");
	fprintf( stderr, "-- Prepare to run...\n");
	fprintf( stderr, "-----------------------------------------------------\n");
	fprintf( stderr, "mdc> #cells = %u = %u x %u x %u\n", (unsigned)( engine.cells.count ), (unsigned)( engine.cells.num.x ),(unsigned)( engine.cells.num.y ),(unsigned)( engine.cells.num.z )  );
	fprintf( stderr, "mdc> cell size= [%g %g %g]\n", engine.cells.length.x, engine.cells.length.y, engine.cells.length.z );
	fprintf( stderr, "mdc> initializing step\n");
	mds_step_init_dt( &engine.step, &engine.frame, engine.dt_sim );
	
	/*==========================================================================
	 *
	 *  run
	 *
	 *========================================================================*/
	if( benchmark )
	{
		double total = 0;
		size_t nbench= strtol( benchmark, NULL, 10);
		if( nbench < 1 ) nbench = 1;
		
		fprintf( stderr, "+++++++++++++++++++++++++++++++++++++++++++++++++\n");
		fprintf( stderr, "++ Benchmarking !!!\n");
		fprintf( stderr, "+++++++++++++++++++++++++++++++++++++++++++++++++\n");
		
		{
			char *filename = mds_string( "bench%u.xyz", (unsigned) np );
			if( filename )
			{
				(void) mds_xyz_save(filename,position,np,filename,0);
				mds_string_delete(filename);
			}
		}
		
		for( i=1; i <= nbench; ++i )
		{
			fprintf( stderr, "[%4u/%4u]  \r", (unsigned)i, (unsigned)nbench); fflush( stderr );
			mds_wtime_start( &chrono );
			mds_dynamics_call( &whole_dynamics, &engine.frame, mds_ljc, &ljc_data);
			total += mds_wtime_query( &chrono );
		}
		fprintf( stderr, "\n");
		total /= nbench;
		total *= 1000;
		fprintf( stderr, "%u => %g ms\n", (unsigned)np, total );
		fprintf( stderr, "V=%g | K=%g | W=%g\n", engine.frame.V, engine.frame.K, engine.frame.W );
		_add_to_benchmark(benchfile, np, total);
		goto MD_END;
		
	}
	
	
	
	/*-------------------------- Initialize Dynamics -------------------------*/
	t_run = t_ini; /** @i_run=0 **/
	fprintf( stderr, "mdc> first whole dynamics...\n");
	mds_wtime_start( &chrono );
	mds_dynamics_call( &whole_dynamics, &engine.frame, mds_ljc, &ljc_data);
	crunch = mds_wtime_query( &chrono );
	fprintf( stderr, "mdc> in %g ms\n", crunch * 1000.0 );
	crunch = 0;
	
	fprintf( stderr, "V=%13.6e | K=%13.6e | W=%13.6e | E=%13.6e\n", 
			engine.frame.V, engine.frame.K, engine.frame.W, engine.frame.K+engine.frame.V); fflush( stderr );
	
	if( !restart )
	{
		(void) save_energy(t_run, &engine.frame, 0);
		if(output)
		{
			if(  mds_file_save( output, &engine.frame, t_run, 0) < 0 )
			{
				mds_perror( stderr, "mds_file_save(%s)", output );
				MD_ERROR;
			}
		}
	}
	
	/*-------------------------- Run !  --------------------------------------*/
	mds_wtime_start( &watch );
	while( t_run < t_max )
	{
		const int save = (0 == ( ++i_run % every ));
		t_run = t_ini + i_run * engine.dt_sim;
		
		mds_wtime_start( &chrono );
		mds_step_predict( & engine.step, &engine.frame);
		mds_dynamics_call( &whole_dynamics, &engine.frame, mds_ljc, &ljc_data);
		mds_step_correct( & engine.step, &engine.frame);
		
		if( rescaleV )
		{
			engine.frame.K = mds_velocity_rescale( velocity, np, (Real)sqrt(K0/engine.frame.K) );
		}
		
		crunch += mds_wtime_query( &chrono );
		
		if( save )
		{
			mds_eta      eta;
			const double fps      = (double)i_run / crunch;         /* "crunch speed", in Frame Per Second */
			const double run_time = mds_wtime_query( &watch );      /* total compute time so far           */
			const double eta_s    = run_time * (t_max/t_run - 1.0); /* Estimated Time of Arrival           */
			mds_eta_compute( &eta, eta_s);
			
			(void) save_energy(t_run, &engine.frame, 1);
			if(output)
			{
				if(  mds_file_save( output, &engine.frame, t_run, 1) < 0 )
				{
					mds_perror( stderr, "mds_file_save(%s)", output );
					MD_ERROR;
				}
			}
			fprintf( stderr, "t= %8.3f |V=%13.6e | K=%13.6e | W=%13.6e | E=%13.6e | <%7.2f> | ETA: %3.0fd%2.0fh%2.0fm%3.1fs \r",
					t_run, engine.frame.V, engine.frame.K, engine.frame.W, engine.frame.K+engine.frame.V, fps, eta.d, eta.h,eta.m, eta.s); 
			fflush( stderr );
		}
		
		
	}
	fprintf( stderr,"\n");
	
	
	/*-------------------------- Done   --------------------------------------*/
	fprintf( stderr, "-----------------------------------------------------\n");
	fprintf( stderr, "-- Done\n");
	fprintf( stderr, "-----------------------------------------------------\n");
	{
		mds_density rho = MDS_DENSITY_INIT;
		if( mds_density_load( &rho, np) < 0 )
		{
			mds_perror(stderr, "mds_density_load" );
			MD_ERROR;
		}
		
		mds_density_build( &rho, position, &engine.frame.box, 2);
		{
			FILE *fp = fopen("density.txt","wb");
			if( fp )
			{
				size_t i;
				for( i=0; i < rho.num_points; ++i )
				{
					fprintf( fp, "%g %g %g\n", rho.X[i], rho.F[i], rho.f[i] );
				}
				fclose(fp);
			}
			
#if 0
			fp = fopen( "density2.txt", "wb" );
			if( fp )
			{
				Real xx = rho.X[0];
				for( ; xx <= rho.X[ rho.num_points - 1]; xx += 0.01 )
				{
					fprintf( fp, "%g %g\n", xx, mds_density_F(&rho,xx) );
				}
				fclose(fp);
			}
#endif		
		}
		mds_density_quit( &rho );
	}
	
	
MD_END:
	MDS_MEMORY();
	mds_options_quit( &options );
	mds_engine_quit( &engine );
	MDS_MEMORY();
	return exit_value;
}
